### Replication code for Intraparty Cleavages and Partisan Attitudes Toward Labor Policy

# This script provides the code required to replicate the article "Intraparty Cleavages and Partisan Attitudes Toward Labor Policy."
# Once the user has downloaded the below raw datasets and added them to the working directory, the entire script can be run from top to bottom.


# NOTE: To run script from top to bottom, user must first download four datasets and add to working directory:

# 1 CCES 2012 (https://dataverse.harvard.edu/dataset.xhtml?persistentId=hdl:1902.1/21447)
# 2 May 2014 BLS Occupational Data (URL to download: https://www.bls.gov/oes/special.requests/oesm14nat.zip)
# 3 Pew March 2015 Survey (http://www.people-press.org/dataset/march-2015-political-survey/) **Log in required**
# 4 General Social Survey 2006 (http://gss.norc.org/get-the-data/stata)

# Packages -----
# Must install the packages below before running script to load them

library(haven)
library(readxl)
library(dplyr)
library(tidyr)
library(ggplot2)
library(purrr)


# Read in raw data ----

cces <- read_dta("CCES12_Common_VV.dta") 
bls <- read_xlsx("oesm14nat/national_M2014_dl.xlsx") 
pew <- read_sav("March15/March15 public.sav")
gss <- read_dta("GSS2006.DTA")

# Correspondence: R objects to tables/figures in paper (NOTE: this section to be used AFTER running full script) -----

# Figure 1:  r_food
# Figure 2:  d_food
# Figure 3:  prop_plot
# Figure 4:  bls_plot
# Figure 5:  ff_pp_plot
# Figure 6:  retail_pp_plot
# Figure 7:  plot_cond_r
# Figure 8:  plot_jobsec_r
# Figure 9:  plot_cond_r_hh
# Figure 10: plot_jobsec_r_hh
# Figure 11: plot_cond_d_hh
# Figure 12: plot_jobsec_d_hh


# To view results from table 1, run the following line of code:
#map(rmods, summary) # must uncomment by removing pound sign to left of code before running

# To view results from table 2, run the following line of code:
#map(dmods, summary) # must uncomment by removing pound sign to left of code before running

# To view results from table 3, run the following line of code:
#map(gssmods, summary) # must uncomment by removing pound sign to left of code before running

# Prepare data: CCES ----

cces12 <- cces %>%
  mutate(
    income5 = ifelse(faminc > 15, NA,
                     ifelse(faminc > 0 & faminc < 3, -2,
                            ifelse(faminc > 2 & faminc < 6, -1,
                                   ifelse(faminc > 5 & faminc < 9, 0,
                                          ifelse(faminc > 8 & faminc < 12, 1,
                                                 ifelse(faminc > 11 & faminc < 16, 2, faminc)))))),
    food_service = ifelse(occupationcat > 14, NA,
                          ifelse(occupationcat == 8, 1, 0)),
    partyid = ifelse(pid3 > 3, NA, pid3),
    rep_id = ifelse(pid7 > 7, NA,
                    ifelse(pid7 > 4 & pid7 < 8, 1, 0)),
    dem_id = ifelse(pid7 > 7, NA,
                    ifelse(pid7 >0 & pid7 < 4, 1, 0))) %>%
  select(income5, food_service, rep_id, dem_id)

# Make income5 a factor variable with middle as reference (treated as set of dummy variables in models, omits middle)
cces12$income5 <- factor(cces12$income5)
levels(cces12$income5) <- c("Lower", "Lower-Middle", "Middle", "Upper-Middle", "Upper")
cces12$income5 <- relevel(cces12$income5, ref = "Middle")


# Prepare data: BLS ----


# Get occ_code, occ_title, h_median 
occ14 <- select(bls, OCC_CODE, OCC_TITLE, H_MEDIAN, A_MEAN)

# Get occupation codes 
occ14 <- occ14 %>%
  filter(OCC_CODE %in% c("33-3050", "33-2011", "25-2000", "35-0000", "41-2010", "51-0000", "53-3021", "53-4000")) %>%
  mutate(
    hr_median_wage = round(as.numeric(H_MEDIAN), 2),
    yr_median_wage = (hr_median_wage*2080)) %>%
  select(occ_code = OCC_CODE, occ_title = OCC_TITLE, hr_median_wage, yr_median_wage, yr_mean = A_MEAN)

# Organize/clean names 
occ14$yr_median_wage <- as.numeric(occ14$yr_median_wage)
occ14[1,2] <- "Teachers"
occ14[4,2] <- "Fast-Food/Food Prep/Service"
occ14[6,2] <- "Manufacturing"
occ14[1,4] <- 55510 # Teachers: only mean available
occ14[9,] <- "Transit"
occ14[10,] <- "Police/Firefighters"
occ14$yr_median_wage <- as.numeric(occ14$yr_median_wage)
occ14[9,4] <- ((occ14[7,4] + occ14[8,4])/2)
occ14[10,4] <- ((occ14[2,4] + occ14[3,4])/2)

df_occ <- occ14 %>%
  filter(occ_title != "Bus Drivers, Transit and Intercity" & occ_title != "Rail Transportation Workers") %>%
  filter(occ_title != "Police Officers" & occ_title != "Firefighters") %>%
  select(occ_title, yr_median_wage)

# Prepare data: Pew ----

dems15 <- pew %>%
  mutate(
    female = ifelse(sex == 2, 1, 0),
    age = ifelse(age > 97, NA, age),
    educ = ifelse(educ2 > 8, NA,
                  ifelse(educ2 >0 & educ2 < 3, 0,
                         ifelse(educ2 > 2 & educ2 < 5, 1,
                                ifelse(educ2 == 5, 2,
                                       ifelse(educ2 > 5 & educ2 < 8, 3,
                                              ifelse(educ2 == 8, 4, educ2)))))),
    latino = ifelse(hisp > 2, NA,
                    ifelse(hisp == 1, 1, 0)),
    black = ifelse(race3m1 > 8, NA,
                   ifelse(race3m1 == 2, 1, 0)),
    income5 = ifelse(income > 9, NA,
                     ifelse(income > 0 & income < 3, -2,
                            ifelse(income > 2 & income  < 5, -1,
                                   ifelse(income > 4 & income < 7, 0,
                                          ifelse(income > 6 & income < 9, 1,
                                                 ifelse(income == 9, 2, income)))))),
    rep_id = ifelse(party > 5, NA,
                    ifelse(party == 1, 1, 0)),
    dem_id = ifelse(party > 5, NA,
                    ifelse(party == 2, 1, 0)),
    ind_id = ifelse(party > 5, NA,
                    ifelse(party == 3, 1, 0)),
    ideology = ifelse(ideo > 5, NA, (max(ideo, na.rm = T) + 1) - ideo),
    union_hh = ifelse(labor2 > 2, NA,
                      ifelse(labor2 == 1, 1, 0)),
    rel_att = ifelse(attend > 6, NA,
                     ifelse(attend > 0 & attend < 3, 1, 0)),
    rep_lean = ifelse(partyln > 2, NA,
                      ifelse(partyln == 1, 1, 0)),
    dem_lean = ifelse(partyln > 2, NA,
                      ifelse(partyln == 2, 1, 0))
  ) %>%
  dplyr::select(female, age, educ, latino, black, income5, rep_id, dem_id, ind_id, ideology, rel_att, union_hh, rep_lean, dem_lean)

atts15 <- pew %>%
  mutate(
    ff = ifelse(q20c > 2, NA,
                ifelse(q20c == 1, 1, 0)),
    retail = ifelse(q20d > 2, NA,
                    ifelse(q20d == 1, 1, 0)),
    pol = ifelse(q20a > 2, NA,
                 ifelse(q20d == 1, 1, 0)),
    teach = ifelse(q20b > 2, NA,
                   ifelse(q20b == 1, 1, 0)),
    factory = ifelse(q20e > 2, NA,
                     ifelse(q20e == 1, 1, 0)),
    transit = ifelse(q20f > 2, NA,
                     ifelse(q20f == 1, 1, 0))
  ) %>%
  select(ff, retail, pol, teach, factory, transit)

# Bind dems and atts
pew15 <- bind_cols(dems15, atts15)

# Make income5 a factor variable with middle as reference
pew15$income5 <- factor(pew15$income5)
levels(pew15$income5) <- c("Lower", "Lower-Middle", "Middle", "Upper-Middle", "Upper")
pew15 <- within(pew15, income5 <- relevel(income5, ref = "Middle"))

# Prepare data: GSS ----

gss06 <- gss %>%
  mutate(
    cond2 = ifelse(unbetter > 5, NA,
                   ifelse(unbetter >=1 & unbetter <= 2, 1, 0)),
    jobsec2 = ifelse(unjobsec > 5, NA,
                     ifelse(unjobsec >=1 & unjobsec <= 2, 1, 0)),
    income5 = ifelse(income06 > 25, NA,
                     ifelse(income06 >= 1 & income06 <= 12, -2,
                            ifelse(income06 > 12 & income06 <= 17, -1,
                                   ifelse(income06 > 17 & income06 <= 20, 0,
                                          ifelse(income06 > 20 & income06 <= 24, 1,
                                                 ifelse(income06 == 25, 2, income06)))))),
    pid_nolean = ifelse(partyid > 6, NA,
                        ifelse(partyid == 0 | partyid == 1, 1,
                               ifelse(partyid > 1 & partyid < 5, 2,
                                      ifelse(partyid == 5 | partyid == 6, 3, partyid)))),
    dem_strong = ifelse(partyid > 6, NA,
                        ifelse(partyid == 0, 1, 0)),
    rep_strong = ifelse(partyid > 6, NA,
                        ifelse(partyid == 6, 1, 0)),
    rep_lean = ifelse(partyid > 6, NA,
                      ifelse(partyid == 4, 1, 0)),
    dem_lean = ifelse(partyid  > 6, NA,
                      ifelse(partyid == 2, 1, 0)),
    educ = ifelse(degree > 4, NA, degree),
    female = ifelse(sex == 2, 1, 0),
    black = ifelse(race == 2, 1, 0),
    latino = ifelse(hispanic > 46, NA,
                    ifelse(hispanic != 1, 1, 0)),
    union_hh = ifelse(union > 4, NA,
                      ifelse(union != 4, 1, 0)),
    ideology = ifelse(polviews > 8, NA,
                      ifelse(polviews == 1, 1,
                             ifelse(polviews == 2 | polviews == 3, 2,
                                    ifelse(polviews == 4 | polviews == 8, 3,
                                           ifelse(polviews == 5 | polviews == 6, 4,
                                                  ifelse(polviews == 7, 5, polviews)))))),
    age = ifelse(age > 100, NA, age),
    rel_att = ifelse(attend > 8, NA,
                     ifelse(attend == 7 | attend == 8, 1, 0))) %>%
  select(cond2, jobsec2, educ, age, female, black, latino, union_hh,
         income5, ideology, rel_att, pid_nolean, dem_strong, rep_strong, rep_lean, dem_lean)

# Make income5 factor for analysis, with reference at middle
gss06$income5 <- factor(gss06$income5)
levels(gss06$income5) <- c("Lower", "Lower-Middle", "Middle", "Upper-Middle", "Upper")
gss06$income5 <- relevel(gss06$income5, ref = "Middle")
# Analyses and Plots: CCES -----

# Republicans
cces12_r <- filter(cces12, rep_id == 1)

# Democrats
cces12_d <- filter(cces12, dem_id == 1)

# Food service
reps_food <- glm(food_service~income5, family = binomial(link = "logit"), data = cces12_r)
dems_food <- glm(food_service~income5, family = binomial(link = "logit"), data = cces12_d)

# Data for predictions
newdata <- data_frame(income5 = c("Lower", "Lower-Middle", "Middle", "Upper-Middle", "Upper"))

# Generate probabilities and CIs: Food service, R & D
pp_r_food <- cbind(newdata, predict(reps_food, newdata, type = "response", se.fit = T))
pp_r_food$upper <- pp_r_food$fit + (1.96*pp_r_food$se.fit)
pp_r_food$lower <- pp_r_food$fit - (1.96*pp_r_food$se.fit)

pp_d_food <- cbind(newdata, predict(dems_food, newdata, type = "response", se.fit = T))
pp_d_food$upper <- pp_d_food$fit + (1.96*pp_d_food$se.fit)
pp_d_food$lower <- pp_d_food$fit - (1.96*pp_d_food$se.fit)

# Reorder factor pp's to correct factor order for plots, R & D
pp_r_food$income5 <- factor(pp_r_food$income5, levels(factor(pp_r_food$income5))[c(1:3, 5, 4)])
pp_d_food$income5 <- factor(pp_d_food$income5, levels(factor(pp_d_food$income5))[c(1:3, 5, 4)])



# Plots
r_food <- ggplot(pp_r_food, aes(x = income5, y = fit)) +
  geom_errorbar(aes(ymax = upper, ymin = lower), width = 0.1) +
  geom_point() +
  labs(title = "Work in Food Preparation Services",
       subtitle = "Republicans",
       y = "Predicted Probability",
       x = "Household Income",
       caption = "Note: Data from CCES 2012.\nResults from logistic regression predicting\nprobability of employment in Food Preparation and Service\nas a function of household income.\nn = 15,675") +
  theme_bw() +
  theme(text = element_text(family = "Times", face = "bold"),
        plot.title = element_text(size = 12),
        axis.text.x = element_text(size = 8),
        plot.caption = element_text(hjust = 0, size = 8)) +
  scale_x_discrete(breaks = c("Lower", "Lower-Middle", "Middle", "Upper-Middle", "Upper"),
                   labels = c("Lower", "Lower-\nMiddle", "Middle", "Upper-\nMiddle", "Upper")) +
  ylim(0, .124)

d_food <- ggplot(pp_d_food, aes(x = income5, y = fit)) +
  geom_errorbar(aes(ymax = upper, ymin = lower), width = 0.1) +
  geom_point() +
  labs(title = "Work in Food Preparation Services",
       subtitle = "Democrats",
       y = "Predicted Probability",
       x = "Household Income",
       caption = "Note: Data from CCES 2012.\nResults from logistic regression predicting\nprobability of employment in Food Preparation and Service\nas a function of household income.\nn = 20,113") +
  theme_bw() +
  theme(text = element_text(family = "Times", face = "bold"),
        plot.title = element_text(size = 12),
        axis.text.x = element_text(size = 8),
        plot.caption = element_text(hjust = 0, size = 8)) +
  scale_x_discrete(breaks = c("Lower", "Lower-Middle", "Middle", "Upper-Middle", "Upper"),
                   labels = c("Lower", "Lower-\nMiddle", "Middle", "Upper-\nMiddle", "Upper")) +
  ylim(0, .124)

# Analyses and Plots: BLS -----

# Caption for BLS plot
cap_note <- paste0(c("\nNote: Data from May 2014 National Occupational Earnings Data, U.S. Bureau of Labor Statistics."),
                   c("\nNational Median Wage Earnings in 2014 = $41,236. Dotted line indicates low-wage threshold, or 2/3 of the median."),
                   c("\nData for police/irefighters is an average of the medians for each occupation."),
                   c("\nData for transit is an avergae of the medians for bus Drivers and rail transit occupations."),
                   c("\nData for teachers refers to mean annual wages."),
                   c("\nBLS occupational codes used: Transit: Mean of Bus Drivers 53-3021, Rail Transport Workers 53-4000;"),
                   c("\nTeachers: 33-2011; Police/Firefighters: Mean of Police Officers 33-3050 and Firefighters 33-2011;"),
                   c("\nManufacturing: 51-0000; Food Prep/Service: 35-0000; Cashiers: 41-2010"))

bls_plot <- ggplot(df_occ) +
  geom_point(aes(x = occ_title, y = yr_median_wage), size = 4) +
  geom_segment(aes(y = 0, x = occ_title, yend = yr_median_wage, xend = occ_title),
               color = "black") +
  geom_hline(yintercept = 27491, linetype = "dotted") +
  labs(title = "Occupations and Wages",
       subtitle = "Left of dotted line: low-wage occupations",
       y = "Annual Median Wage Earnings ($)",
       x = NULL,
       caption = cap_note) +
  scale_y_continuous(breaks = c(0, 10000, 20000, 30000, 40000, 50000)) +
  theme_bw() +
  theme(text = element_text(family = "Times", face = "bold")) +
  theme(plot.caption = element_text(hjust = 0, size = 8)) +
  coord_flip()

# Analyses and Plots: Pew -----

# Proportions support whole Pew Sample
attsplot15 <- gather(atts15, occ, support, 1:6)

# Means and SEs
props <- attsplot15 %>%
  filter(!is.na(occ)) %>%
  filter(!is.na(support)) %>%
  group_by(occ) %>%
  summarize(avg= mean(support, na.rm = T),
            N = n(),
            sd = sd(support, na.rm = T),
            se = sd/sqrt(N))

# Names
props <- mutate(props, occ_names = ifelse(occ == "ff", "Fast-Food",
                                          ifelse(occ == "pol", "Police/Firefighter",
                                                 ifelse(occ == "teach", "Teacher",
                                                        ifelse(occ == "retail", "Supermarket/Retail",
                                                               ifelse(occ == "transit", "Transit",
                                                                      ifelse(occ == "factory", "Factory/Manufact.", occ)))))))
# Plot
prop_plot <- ggplot(props, aes(x = occ_names, y = avg)) +
  geom_col(fill = "gray") + theme_bw() +
  geom_errorbar(aes(ymin = avg-se, ymax = avg+se), width = 0.2) +
  labs(title = "Proportion of Support by Sector Unionization",
       subtitle = "Full Sample",
       x = NULL,
       y = "Percentage Support",
       caption = "Views on unionization in different industries. Note: Pew Political Survey 2015") +
  theme(text = element_text(family = "Times")) + 
  scale_y_continuous(labels = scales::percent) + coord_flip()

# Models for Pew results section

# Filter for only Republicans (& leaners)
pew15_r <- filter(pew15, rep_id == 1 | rep_lean == 1)

# Republican models
pol_mod_r <- glm(pol ~ income5  + union_hh + age + female + black + latino +
                   educ + ideology + rel_att + rep_id,
                 family = binomial("logit"), data = pew15_r)
teach_mod_r <- glm(teach ~ income5  + union_hh + age + female + black + latino +
                   educ + ideology + rel_att + rep_id,
                 family = binomial("logit"), data = pew15_r)
ff_mod_r <- glm(ff ~ income5  + union_hh + age + female + black + latino +
                   educ + ideology + rel_att + rep_id,
                 family = binomial("logit"), data = pew15_r)
retail_mod_r <- glm(retail ~ income5  + union_hh + age + female + black + latino +
                   educ + ideology + rel_att + rep_id,
                 family = binomial("logit"), data = pew15_r)
factor_mod_r <- glm(factory ~ income5  + union_hh + age + female + black + latino +
                   educ + ideology + rel_att + rep_id,
                 family = binomial("logit"), data = pew15_r)
transit_mod_r <- glm(transit ~ income5  + union_hh + age + female + black + latino +
                   educ + ideology + rel_att + rep_id,
                 family = binomial("logit"), data = pew15_r)

rmods <- list(ff_mod_r, retail_mod_r, teach_mod_r, pol_mod_r, factor_mod_r, transit_mod_r)

# Filter for only Democrats (& leaners)
pew15_d <- filter(pew15, dem_id == 1 | dem_lean == 1)

# Democratic models
pol_mod_d <- glm(pol ~ income5  + union_hh + age + female + black + latino +
                   educ + ideology + rel_att + dem_id,
                 family = binomial("logit"), data = pew15_d)
teach_mod_d <- glm(teach ~ income5  + union_hh + age + female + black + latino +
                     educ + ideology + rel_att + dem_id,
                   family = binomial("logit"), data = pew15_d)
ff_mod_d <- glm(ff ~ income5  + union_hh + age + female + black + latino +
                  educ + ideology + rel_att + dem_id,
                family = binomial("logit"), data = pew15_d)
retail_mod_d <- glm(retail ~ income5  + union_hh + age + female + black + latino +
                      educ + ideology + rel_att + dem_id,
                    family = binomial("logit"), data = pew15_d)
factor_mod_d <- glm(factory ~ income5  + union_hh + age + female + black + latino +
                      educ + ideology + rel_att + dem_id,
                    family = binomial("logit"), data = pew15_d)
transit_mod_d <- glm(transit ~ income5  + union_hh + age + female + black + latino +
                       educ + ideology + rel_att + dem_id,
                     family = binomial("logit"), data = pew15_d)

dmods <- list(ff_mod_d, retail_mod_d, teach_mod_d, pol_mod_d, factor_mod_d, transit_mod_d)


# Predicted probabilities for Republicans by income

# Generate new data for predictions
newframe <- data_frame(income5 = c("Lower", "Lower-Middle", "Middle", "Upper-Middle", "Upper"),
                       age = median(pew15_r$age, na.rm = T),
                       educ = median(pew15_r$educ, na.rm = T),
                       black = median(pew15_r$black, na.rm = T),
                       latino = median(pew15_r$latino, na.rm = T),
                       female = median(pew15_r$female, na.rm = T),
                       union_hh = median(pew15_r$union_hh, na.rm = T),
                       ideology = median(pew15_r$ideology, na.rm = T),
                       rep_id = median(pew15_r$rep_id, na.rm = T),
                       rel_att = median(pew15_r$rel_att, na.rm = T))

# Predicted probabilities
newframe_ff <- cbind(newframe, predict(ff_mod_r, newframe, type = "response", se.fit = T))
newframe_retail <- cbind(newframe, predict(retail_mod_r, newframe, type = "response", se.fit = T))

# CIs
newframe_ff$lower <- newframe_ff$fit - (1.96*newframe_ff$se.fit)
newframe_ff$upper <- newframe_ff$fit + (1.96*newframe_ff$se.fit)

newframe_retail$lower <- newframe_retail$fit - (1.96*newframe_retail$se.fit)
newframe_retail$upper <- newframe_retail$fit + (1.96*newframe_retail$se.fit)

# Correct factor for plot
newframe_ff$income5 <- factor(newframe_ff$income5, levels(factor(newframe_ff$income5))[c(1, 2, 3, 5, 4)])
newframe_retail$income5 <- factor(newframe_retail$income5, levels(factor(newframe_retail$income5))[c(1, 2, 3, 5, 4)])

# FF and retail plots
ff_pp_plot <- ggplot(newframe_ff, aes(x = income5, y = fit)) +
  geom_errorbar(aes(ymax = upper, ymin = lower), width = 0.1) +
  geom_point() +
  labs(title = "Fast-Food Worker Unionization",
       subtitle = "Republican Support",
       y = "Predicted Probability",
       x = "Household Income") + theme_bw() +
  theme(text = element_text(family = "Times", face = "bold"),
        plot.title = element_text(size = 12),
        axis.text.x = element_text(size = 8)) +
  scale_x_discrete(breaks = c("Lower", "Lower-Middle", "Middle", "Upper-Middle", "Upper"),
                   labels = c("Lower", "Lower-\nMiddle", "Middle", "Upper-\nMiddle", "Upper")) 

retail_pp_plot <- ggplot(newframe_retail, aes(x = income5, y = fit)) +
  geom_errorbar(aes(ymax = upper, ymin = lower), width = 0.1) +
  geom_point() +
  labs(title = "Retail Worker Unionization",
       subtitle = "Republican Support",
       y = "Predicted Probability",
       x = "Household Income") + theme_bw() +
  theme(text = element_text(family = "Times", face = "bold"),
        plot.title = element_text(size = 12),
        axis.text.x = element_text(size = 8)) +
  scale_x_discrete(breaks = c("Lower", "Lower-Middle", "Middle", "Upper-Middle", "Upper"),
                   labels = c("Lower", "Lower-\nMiddle", "Middle", "Upper-\nMiddle", "Upper"))

# Calculations for police/firefighter and teachers, results in text not plots
newframe_pol <- cbind(newframe, predict(pol_mod_r, newframe, type = "response", se.fit = T))
newframe_pol$lower <- newframe_pol$fit - (1.96*newframe_pol$se.fit)
newframe_pol$upper <- newframe_pol$fit + (1.96*newframe_pol$se.fit)

newframe_teach <- cbind(newframe, predict(teach_mod_r, newframe, type = "response", se.fit = T))
newframe_teach$lower <- newframe_teach$fit - (1.96*newframe_teach$se.fit)
newframe_teach$upper <- newframe_teach$fit + (1.96*newframe_teach$se.fit)

# Predicted probabilities for Republican union hhs
# New data to generate predictions
newframe_hh <- data_frame(income5 = c("Middle"),
                          age = median(pew15_r$age, na.rm = T),
                          educ = median(pew15_r$educ, na.rm = T),
                          black = median(pew15_r$black, na.rm = T),
                          latino = median(pew15_r$latino, na.rm = T),
                          female = median(pew15_r$female, na.rm = T),
                          union_hh = c(0, 1),
                          ideology = median(pew15_r$ideology, na.rm = T),
                          rep_id = median(pew15$rep_id, na.rm = T),
                          rel_att = median(pew15_r$rel_att, na.rm = T))

## FF and retail, union hhs
newframe_ff <- cbind(newframe_hh, predict(ff_mod_r, newframe_hh, type = "response", se.fit = T))
newframe_ff$lower <- newframe_ff$fit - (1.96*newframe_ff$se.fit)
newframe_ff$upper <- newframe_ff$fit + (1.96*newframe_ff$se.fit)

newframe_retail <- cbind(newframe_hh, predict(retail_mod_r, newframe_hh, type = "response", se.fit = T))
newframe_retail$lower <- newframe_retail$fit - (1.96*newframe_retail$se.fit)
newframe_retail$upper <- newframe_retail$fit + (1.96*newframe_retail$se.fit)


# Analyses and Plots: GSS -----

# Republicans
rep06 <- filter(gss06, pid_nolean == 3 | rep_lean == 1)

cond_rep <- glm(cond2 ~ income5  + union_hh + age + female + black + latino +
                  educ + ideology + rel_att + rep_strong,
                family = binomial("logit"), data = rep06)
jobsec_rep <- glm(jobsec2 ~ income5  + union_hh + age + female + black + latino +
                    educ + ideology + rel_att + rep_strong,
                  family = binomial("logit"), data = rep06)

# Democrats
dem06 <- filter(gss06, pid_nolean == 1 | dem_lean == 1)

cond_dem <- glm(cond2 ~ income5  + union_hh + age + female + black + latino +
                  educ + ideology + rel_att + dem_strong,
                family = binomial("logit"), data = dem06)
jobsec_dem <- glm(jobsec2 ~ income5  + union_hh + age + female + black + latino +
                    educ + ideology + rel_att + dem_strong,
                  family = binomial("logit"), data = dem06)

gssmods <- list(cond_rep, jobsec_rep, cond_dem, jobsec_dem)

# Generate newdata for predictions, Republicans, income
newdata_r <- data_frame(income5 = c("Lower", "Lower-Middle", "Middle", "Upper-Middle", "Upper"),
                        age = median(rep06$age, na.rm = T),
                        educ = median(rep06$educ, na.rm = T),
                        black = median(rep06$black, na.rm = T),
                        latino = median(rep06$latino, na.rm = T),
                        female = median(rep06$female, na.rm = T),
                        union_hh = median(rep06$union_hh, na.rm = T),
                        ideology = median(rep06$ideology, na.rm = T),
                        rep_strong = median(rep06$rep_strong, na.rm = T),
                        rel_att = median(rep06$rel_att, na.rm = T))

# Predicted probabilities & CIs
pp_cond_r <- cbind(newdata_r, predict(cond_rep, newdata_r, type = "response", se.fit = T))
pp_cond_r$lower <- pp_cond_r$fit - (1.96*pp_cond_r$se.fit)
pp_cond_r$upper <- pp_cond_r$fit + (1.96*pp_cond_r$se.fit)

pp_jobsec_r <- cbind(newdata_r, predict(jobsec_rep, newdata_r, type = "response", se.fit = T))
pp_jobsec_r$lower <- pp_jobsec_r$fit - (1.96*pp_jobsec_r$se.fit)
pp_jobsec_r$upper <- pp_jobsec_r$fit + (1.96*pp_jobsec_r$se.fit)

# Income as factor ordered for plot
pp_cond_r$income5 <- factor(pp_cond_r$income5, levels(factor(pp_cond_r$income5))[c(1, 2, 3, 5, 4)])
pp_jobsec_r$income5 <- factor(pp_jobsec_r$income5, levels(factor(pp_jobsec_r$income5))[c(1, 2, 3, 5, 4)])

# Plots, Republicans by income
plot_cond_r <- ggplot(pp_cond_r, aes(x = income5, y = fit)) + geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.1) +
  labs(title = "Unions Improve Working Conditions",
       subtitle = "Republican Agreement",
       y = "Predicted Probability",
       x = "Household Income") +
  theme_bw() +
  theme(text = element_text(family = "Times", face = "bold"),
        plot.title = element_text(size = 12),
        axis.text.x = element_text(size = 8)) +
  scale_x_discrete(breaks = c("Lower", "Lower-Middle", "Middle", "Upper-Middle", "Upper"),
                   labels = c("Lower", "Lower-\nMiddle", "Middle", "Upper-\nMiddle", "Upper")) 

plot_jobsec_r <- ggplot(pp_jobsec_r, aes(x = income5, y = fit)) + geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.1) +
  labs(title = "Unions Improve Worker Job Security",
       subtitle = "Republican Agreement",
       y = "Predicted Probability",
       x = "Household Income") +
  theme_bw() +
  theme(text = element_text(family = "Times", face = "bold"),
        plot.title = element_text(size = 12),
        axis.text.x = element_text(size = 8)) +
  scale_x_discrete(breaks = c("Lower", "Lower-Middle", "Middle", "Upper-Middle", "Upper"),
                   labels = c("Lower", "Lower-\nMiddle", "Middle", "Upper-\nMiddle", "Upper")) 

# Generate data for union hhs predicted probabilities
# union hh, R
newdata_r_hh <- data_frame(income5 = c("Middle"),
                           age = median(rep06$age, na.rm = T),
                           educ = median(rep06$educ, na.rm = T),
                           black = median(rep06$black, na.rm = T),
                           latino = median(rep06$latino, na.rm = T),
                           female = median(rep06$female, na.rm = T),
                           union_hh = c(0, 1),
                           ideology = median(rep06$ideology, na.rm = T),
                           rep_strong = median(rep06$rep_strong, na.rm = T),
                           rel_att = median(rep06$rel_att, na.rm = T))

# union hh, D
newdata_d_hh <- data_frame(income5 = c("Middle"),
                           age = median(dem06$age, na.rm = T),
                           educ = median(dem06$educ, na.rm = T),
                           black = median(dem06$black, na.rm = T),
                           latino = median(dem06$latino, na.rm = T),
                           female = median(dem06$female, na.rm = T),
                           union_hh = c(0, 1),
                           ideology = median(dem06$ideology, na.rm = T),
                           dem_strong = median(dem06$dem_strong, na.rm = T),
                           rel_att = median(dem06$rel_att, na.rm = T))

# Predicted probabilities, union hh, Republicans
pp_cond_r_hh <- cbind(newdata_r_hh, predict(cond_rep, newdata_r_hh, type = "response", se.fit = T))
pp_cond_r_hh$lower <- pp_cond_r_hh$fit - (1.96*pp_cond_r_hh$se.fit)
pp_cond_r_hh$upper <- pp_cond_r_hh$fit + (1.96*pp_cond_r_hh$se.fit)

pp_jobsec_r_hh <- cbind(newdata_r_hh, predict(jobsec_rep, newdata_r_hh, type = "response", se.fit = T))
pp_jobsec_r_hh$lower <- pp_jobsec_r_hh$fit - (1.96*pp_jobsec_r_hh$se.fit)
pp_jobsec_r_hh$upper <- pp_jobsec_r_hh$fit + (1.96*pp_jobsec_r_hh$se.fit)

# Predicted probabilities, union hh, Democrats
pp_cond_d_hh <- cbind(newdata_d_hh, predict(cond_dem, newdata_d_hh, type = "response", se.fit = T))
pp_cond_d_hh$lower <- pp_cond_d_hh$fit - (1.96*pp_cond_d_hh$se.fit)
pp_cond_d_hh$upper <- pp_cond_d_hh$fit + (1.96*pp_cond_d_hh$se.fit)

pp_jobsec_d_hh <- cbind(newdata_d_hh, predict(jobsec_dem, newdata_d_hh, type = "response", se.fit = T))
pp_jobsec_d_hh$lower <- pp_jobsec_d_hh$fit - (1.96*pp_jobsec_d_hh$se.fit)
pp_jobsec_d_hh$upper <- pp_jobsec_d_hh$fit + (1.96*pp_jobsec_d_hh$se.fit)

# Predicted probabilities: views on unions and working conditions for lower-income, extremely conservative Republican
newdata_wc_r <- data_frame(income5 = c("Lower"),
                        age = median(rep06$age, na.rm = T),
                        educ = median(rep06$educ, na.rm = T),
                        black = median(rep06$black, na.rm = T),
                        latino = median(rep06$latino, na.rm = T),
                        female = median(rep06$female, na.rm = T),
                        union_hh = median(rep06$union_hh, na.rm = T),
                        ideology = 5,
                        rep_strong = median(rep06$rep_strong, na.rm = T),
                        rel_att = median(rep06$rel_att, na.rm = T))

# Predicted probabilities: views on unions and working conditions for upper-middle-income, extremely liberal Democrat
newdata_wc_d <- data_frame(income5 = c("Upper-Middle"),
                        age = median(dem06$age, na.rm = T),
                        educ = median(dem06$educ, na.rm = T),
                        black = median(dem06$black, na.rm = T),
                        latino = median(dem06$latino, na.rm = T),
                        female = median(dem06$female, na.rm = T),
                        union_hh = median(dem06$union_hh, na.rm = T),
                        ideology = 1,
                        dem_strong = median(dem06$dem_strong, na.rm = T),
                        rel_att = median(dem06$rel_att, na.rm = T))


pp_con_r <- cbind(newdata_wc_r, predict(cond_rep, newdata_wc_r, type = "response", se.fit = T))
pp_con_r$lower <- pp_con_r$fit - (1.96*pp_con_r$se.fit)
pp_con_r$upper <- pp_con_r$fit + (1.96*pp_con_r$se.fit)

pp_lib_d <- cbind(newdata_wc_d, predict(cond_dem, newdata_wc_d, type = "response", se.fit = T))
pp_lib_d$lower <- pp_lib_d$fit - (1.96*pp_lib_d$se.fit)
pp_lib_d$upper <- pp_lib_d$fit + (1.96*pp_lib_d$se.fit)

# Remove pound sign to uncomment and print results
#pp_con_r
#pp_lib_d

## Plots, Republicans, union hh
plot_cond_r_hh <- ggplot(pp_cond_r_hh, aes(x = union_hh, y = fit)) + geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.1) +
  scale_x_continuous(breaks = c(0, 1), labels = c("Non-Union", "Union")) +
  labs(title = "Unions Improve Working Conditions",
       subtitle = "Republican Agreement",
       y = "Predicted Probability",
       x = "Household") +
  theme_bw() +
  theme(text = element_text(family = "Times", face = "bold"),
        plot.title = element_text(size = 12),
        axis.text.x = element_text(size = 8)) 

plot_jobsec_r_hh <- ggplot(pp_jobsec_r_hh, aes(x = union_hh, y = fit)) + geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.1) +
  scale_x_continuous(breaks = c(0, 1), labels = c("Non-Union", "Union")) +
  labs(title = "Unions Improve Job Security",
       subtitle = "Republican Agreement",
       y = "Predicted Probability",
       x = "Household") +
  theme_bw() +
  theme(text = element_text(family = "Times", face = "bold"),
        plot.title = element_text(size = 12),
        axis.text.x = element_text(size = 8)) 

## Plots, Democrats, union hhs
plot_cond_d_hh <- ggplot(pp_cond_d_hh, aes(x = union_hh, y = fit)) + geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.1) +
  scale_x_continuous(breaks = c(0, 1), labels = c("Non-Union", "Union")) +
  labs(title = "Unions Improve Working Conditions",
       subtitle = "Democratic Agreement",
       y = "Predicted Probability",
       x = "Household") +
  theme_bw() +
  theme(text = element_text(family = "Times", face = "bold"),
        plot.title = element_text(size = 12),
        axis.text.x = element_text(size = 8)) 

plot_jobsec_d_hh <- ggplot(pp_jobsec_d_hh, aes(x = union_hh, y = fit)) + geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.1) +
  scale_x_continuous(breaks = c(0, 1), labels = c("Non-Union", "Union")) +
  labs(title = "Unions Improve Job Security",
       subtitle = "Democratic Agreement",
       y = "Predicted Probability",
       x = "Household") +
  theme_bw() +
  theme(text = element_text(family = "Times", face = "bold"),
        plot.title = element_text(size = 12),
        axis.text.x = element_text(size = 8)) 


# Tables in paper (must uncomment and be used with knitr and texreg package) ----

# Labels for Pew table (Tables 1 & 2)

mod_names <- c("Fast-Food",
               "Retail",
               "Teachers",
               "Police/Firefighters",
               "Factory/Manufact.",
               "Transit")
ivs_d <- c("Intercept",
           "Lower Income",
           "Lower-Middle Income",
           "Upper-Middle Income",
           "Upper Income",
           "Union Household",
           "Age",
           "Female",
           "Black",
           "Latino",
           "Education",
           "Ideology",
           "Religious Attendance", 
           "Strong Democrat")
ivs_r <- c("Intercept",
           "Lower Income",
           "Lower-Middle Income",
           "Upper-Middle Income",
           "Upper Income",
           "Union Household",
           "Age",
           "Female",
           "Black",
           "Latino",
           "Education",
           "Ideology",
           "Religious Attendance", 
           "Strong Republican")

# Labels for GSS table (Table 3)

gssmod_names <- c("Reps: Work Conditions",
                  "Reps: Job Security",
                  "Dems: Work Conditions",
                  "Dems: Job Security")

ivs_gss <- c("Intercept",
             "Lower Income",
             "Lower-Middle Income",
             "Upper-Middle Income",
             "Upper Income",
             "Union Household",
             "Age",
             "Female",
             "Black",
             "Latino",
             "Education", 
             "Ideology",
             "Religious Attendance", 
             "Strong Republican",
             "Strong Democrat")


## Table 1: Republican Partisan Attitudes Toward Unionization

# texreg(rmods, digits = 2, caption = "Republican Partisan Attitudes Toward Unionization",
#        center = T, scalebox = .85, caption.above = T,
#        custom.coef.names = ivs_r,
#        reorder.coef = c(6, 2:5, 7:14, 1),
#        custom.model.names = mod_names)

## Table 2: Democratic Partisan Attitudes Toward Unionization

# texreg(dmods, digits = 2, caption = "Democratic Partisan Attitudes Toward Unionization",
#        center = T, caption.above = T, scalebox = .85,
#        custom.coef.names = ivs_d,
#        reorder.coef = c(6, 2:5, 7:14, 1),
#        custom.model.names = mod_names)

## Table 3: Partisan Attitudes Toward Unions in the Workplace

# texreg(gssmods, digits = 2, caption = "Partisan Attitudes Toward Unions in the Workplace",
#        center = T, caption.above = T, scalebox = .85,
#        custom.coef.names = ivs_gss,
#        reorder.coef = c(6, 2:5, 7:15, 1),
#        custom.model.names = gssmod_names)
